#################################
# Representative Matching  
#################################

rm(list=ls())

# load packages
library(devtools)
library(slam)
library(designmatch) 
library(foreign) 
library(gurobi)

# load data
load("anticorr_replication_data.RData")

###########################
#### Prepare variables #### 
###########################

# age
d$age_resp = d$Q1
table(d$age_resp, exclude = NULL)

# female
d$female_resp = 0
d$female_resp[d$Q2=="Mujer"] = 1
table(d$female_resp, exclude = NULL)

# education 
d$more_high_school = 0
d$more_high_school[d$a1=="13 años: 1er año de educación terciaria o superior universitaria"] = 1
d$more_high_school[d$a1=="14 años: 2do año de educación terciaria o superior universitaria"] = 1
d$more_high_school[d$a1=="15 años: 3er año de educación terciaria o superior universitaria"] = 1
d$more_high_school[d$a1=="16 años: 4to año de educación terciaria o superior universitaria"] = 1
d$more_high_school[d$a1=="17 años: 5to año de educación terciaria o superior universitaria"] = 1
d$more_high_school[d$a1=="18 años o más: 6to año de educación terciaria o superior universitaria"] = 1
table(d$more_high_school)

##############################################
#### Matching concrete anti-corruption #######
##############################################

# treatment concrete
d$concrete_anti = 0
d$concrete_anti[d$atributo4=="•\tSu partido propone combatir la corrupción, y ofrece medidas concretas anticorrupción como transparentar la gestión pública."] = 1
d$concrete_anti[d$atributo4=="•\tSu partido propone combatir la corrupción, pero no ofrece medidas concretas anticorrupción."] = 2
table(d$concrete_anti, exclude = NULL)
d <- d[d$concrete_anti != 2,]

# sort by treatment
d = d[order(d$concrete_anti, decreasing=TRUE), ]

# treatment indicator
t_ind = d$concrete_anti

# moment balance
mom_covs = cbind(d$age_resp,d$female_resp,d$more_high_school) # age, sex, education
mom_targets = c(28, 0.5, 0.25) # target distribution
mom_tols = absstddif(mom_covs, t_ind, .05)
mom = list(covs = mom_covs, tols = mom_tols, targets = mom_targets)

# solver paramaters
t_max = 60*30
name = "gurobi"
approximate = 0
solver = list(name = name, t_max = t_max, approximate = approximate, round_cplex = 0, trace = 1)

# find optimal matches
out_1 = cardmatch(t_ind, mom = mom, solver = solver)

# numbers
t_id_1 = out_1$t_id
c_id_1 = out_1$c_id
length(t_id_1) # [1] 5133
length(c_id_1) # [1] 5133

# save matched sample
d_match_concrete = d[c(out_1$t_id, out_1$c_id), ] 
save(d_match_concrete, file = "d_match_concrete.Rdata")

########################################################
#### Replicate concrete anti-corruption w/matched sample 
########################################################

# make concrete a factor
table(d_match_concrete$atributo4)
d_match_concrete$concrete_anti = 0
d_match_concrete$concrete_anti[d_match_concrete$atributo4=="•\tSu partido propone combatir la corrupción, y ofrece medidas concretas anticorrupción como transparentar la gestión pública."] = 1
table(d_match_concrete$concrete_anti, exclude = NULL)
d_match_concrete$concrete_anti = factor(d_match_concrete$concrete_anti, 0:1, c("No Anticorr.", "Concrete Anticorr."))

# analysis 1: vote - concrete
m1 <- outcome_q1 ~ concrete_anti + gender + ideology + party + bribery + embezzlement + competence + age + education
amces1 <- cj(d_match_concrete, m1, id = ~id)
head(amces1[c("feature", "level", "estimate", "std.error", "p", "lower", "upper")], 20L)

# analysis 2: perceived effectiveness - concrete
m2 <- outcome_q4 ~ concrete_anti + gender + ideology + party + bribery + embezzlement + competence + age + education
amces2 <- cj(d_match_concrete, m2, id = ~id)
head(amces2[c("feature", "level", "estimate", "std.error", "p", "lower", "upper")], 20L)

#############################################
#### Matching vague anti-corruption #########
#############################################

# treatment vague
table(d$atributo4)
d$vague_anti = 0
d$vague_anti[d$atributo4=="•\tSu partido propone combatir la corrupción, pero no ofrece medidas concretas anticorrupción."] = 1
d$vague_anti[d$atributo4=="•\tSu partido propone combatir la corrupción, y ofrece medidas concretas anticorrupción como transparentar la gestión pública."] = 2
table(d$vague_anti, exclude = NULL)
d <- d[d$vague_anti != 2,]

# sort by treatment
d = d[order(d$vague_anti, decreasing=TRUE), ]

# treatment indicator
t_ind = d$vague_anti

# moment balance
mom_covs = cbind(d$age_resp,d$female_resp,d$more_high_school) # age, sex, education
mom_targets = c(28, 0.5, 0.25) # target distribution
mom_tols = absstddif(mom_covs, t_ind, .05) # constrain at most .05 standard deviations apart
mom = list(covs = mom_covs, tols = mom_tols, targets = mom_targets)

# solver paramaters
t_max = 60*30
name = "gurobi"
approximate = 0
solver = list(name = name, t_max = t_max, approximate = approximate, round_cplex = 0, trace = 1)

# find optimal matches
out_1 = cardmatch(t_ind, mom = mom, solver = solver)

# numbers
t_id_1 = out_1$t_id
c_id_1 = out_1$c_id
length(t_id_1) # [1] 3164
length(c_id_1) # [1] 3164

# save matched sample
d_match_vague = d[c(out_1$t_id, out_1$c_id), ] 
save(d_match_vague, file = "d_match_vague.Rdata")

#####################################################
### Replicate vague anti-corruption w/matched sample 
#####################################################

# make vague a factor
d_match_vague$vague_anti = 0
d_match_vague$vague_anti[d_match_vague$atributo4=="•\tSu partido propone combatir la corrupción, pero no ofrece medidas concretas anticorrupción."] = 1
d_match_vague$vague_anti[d_match_vague$atributo4=="•\tSu partido propone combatir la corrupción, y ofrece medidas concretas anticorrupción como transparentar la gestión pública."] = 2
table(d_match_vague$vague_anti, exclude = NULL)
d_match_vague$vague_anti = factor(d_match_vague$vague_anti, 0:1, c("No Anticorr.", "Vague Anticorr."))

# analysis 1: vote - vague
m1 <- outcome_q1 ~ vague_anti + gender + ideology + party + bribery + embezzlement + competence + age + education
amces1 <- cj(d_match_vague, m1, id = ~id)
head(amces1[c("feature", "level", "estimate", "std.error", "p", "lower", "upper")], 20L)

# analysis 2: perceived effectiveness - vague
m2 <- outcome_q4 ~ vague_anti + gender + ideology + party + bribery + embezzlement + competence + age + education
amces2 <- cj(d_match_vague, m2, id = ~id)
head(amces2[c("feature", "level", "estimate", "std.error", "p", "lower", "upper")], 20L)

##########################
#### Matching bribery ####
##########################

# treatment
table(d$atributo1)
d$bribery = 0
d$bribery[d$atributo1=="•\tSi ha recibido sobornos a cambio de contratos de obra pública."] = 1
table(d$bribery, exclude = NULL)

# sort by treatment
d = d[order(d$bribery, decreasing=TRUE), ]

# treatment indicator
t_ind = d$bribery

# moment balance
mom_covs = cbind(d$age_resp,d$female_resp,d$more_high_school) # age, sex, education
mom_targets = c(28, 0.5, 0.25) # target distribution
mom_tols = absstddif(mom_covs, t_ind, .05)
mom = list(covs = mom_covs, tols = mom_tols, targets = mom_targets)

# solver paramaters
t_max = 60*30
name = "gurobi"
approximate = 0
solver = list(name = name, t_max = t_max, approximate = approximate, round_cplex = 0, trace = 1)

# find optimal matches
out_1 = cardmatch(t_ind, mom = mom, solver = solver)

# numbers
t_id_1 = out_1$t_id
c_id_1 = out_1$c_id
length(t_id_1) # [1] 6826
length(c_id_1) # [1] 6826

# save matched sample
d_match_bribery = d[c(out_1$t_id, out_1$c_id), ] 
save(d_match_bribery, file = "d_match_bribery.RData")

############################
# Replicate bribery analysis 
############################

# make bribery a factor
d_match_bribery$bribery = NA
d_match_bribery$bribery[d_match_bribery$atributo1=="•\tNo ha recibido sobornos a cambio de contratos de obra pública."] = 0
d_match_bribery$bribery[d_match_bribery$atributo1=="•\tSi ha recibido sobornos a cambio de contratos de obra pública."] = 1
table(d_match_bribery$bribery, exclude = NULL)
d_match_bribery$bribery = factor(d_match_bribery$bribery, 0:1, c("No Bribes", "Bribes"))

# analysis 3: vote by bribery
amces3 <- cj(d_match_bribery, outcome_q1 ~ anticorruption + gender + ideology + party + embezzlement + competence + age + education, id = ~id, estimate = "amce", by = ~bribery)
diff_amces3 <- cj(d_match_bribery, outcome_q1 ~ anticorruption + gender + ideology + party + embezzlement + competence + age + education, id = ~id, estimate = "amce_diff", by = ~bribery)
plot3 = plot(rbind(amces3, diff_amces3)) + ggplot2::facet_wrap(~BY, ncol = 3L) + ggtitle("Change in Electoral Support") + theme(plot.title = element_text(hjust = 0.5))
plot3

# figure A2
plot3 <- plot3 + 
  scale_color_manual(values = rep("black", 9))+
  theme(legend.position = "none")
plot3

# analysis 4: perceived effectiveness by bribery
amces4 <- cj(d_match_bribery, outcome_q4 ~ anticorruption + gender + ideology + party + embezzlement + competence + age + education, id = ~id, estimate = "amce", by = ~bribery)
diff_amces4 <- cj(d_match_bribery, outcome_q4 ~ anticorruption + gender + ideology + party + embezzlement + competence + age + education, id = ~id, estimate = "amce_diff", by = ~bribery)
plot4 = plot(rbind(amces4, diff_amces4)) + ggplot2::facet_wrap(~BY, ncol = 3L) + ggtitle("Change in Perceived Effectiveness") + theme(plot.title = element_text(hjust = 0.5))

# figure A1
plot4 <- plot4 + 
  scale_color_manual(values = rep("black", 9))+
  theme(legend.position = "none")
plot4

#########################
#### Matching gender ####
#########################

# treatment
d$gender = 0
d$gender[d$atributo7=="•\tMujer"] = 1
table(d$gender, exclude = NULL)

# sort by treatment
d = d[order(d$gender, decreasing=TRUE), ]

# treatment indicator
t_ind = d$gender

# moment balance
mom_covs = cbind(d$age_resp,d$female_resp,d$more_high_school) # age, sex, education
mom_targets = c(28, 0.5, 0.25) # target distribution
mom_tols = absstddif(mom_covs, t_ind, .05)
mom = list(covs = mom_covs, tols = mom_tols, targets = mom_targets)

# solver paramaters
t_max = 60*30
name = "gurobi"
approximate = 0
solver = list(name = name, t_max = t_max, approximate = approximate, round_cplex = 0, trace = 1)

# find optimal matches
out_1 = cardmatch(t_ind, mom = mom, solver = solver)

# numbers
t_id_1 = out_1$t_id
c_id_1 = out_1$c_id
length(t_id_1) # [1] 6926
length(c_id_1) # [1] 6926

# save matched sample
d_match_gender = d[c(out_1$t_id, out_1$c_id), ] 
save(d_match_gender, file = "d_match_gender.RData")

############################## 
#### Replicate gender analysis
##############################

# make gender a factor
d_match_gender$gender = NA
d_match_gender$gender[d_match_gender$atributo7=="•\tHombre"] = 0
d_match_gender$gender[d_match_gender$atributo7=="•\tMujer"] = 1
table(d_match_gender$gender, exclude = NULL)
d_match_gender$gender = factor(d_match_gender$gender, 0:1, c("Man", "Woman"))

# analysis 5: perceived effectiveness by gender
amces5 <- cj(d_match_gender, outcome_q4 ~ anticorruption + ideology + party + bribery + embezzlement + competence + age + education, id = ~id, estimate = "amce", by = ~gender)
diff_amces5 <- cj(d_match_gender, outcome_q4 ~ anticorruption + ideology + party + bribery + embezzlement + competence + age + education, id = ~id, estimate = "amce_diff", by = ~gender)
plot5 = plot(rbind(amces5, diff_amces5)) + ggplot2::facet_wrap(~BY, ncol = 3L) + ggtitle("Change in Perceived Effectiveness: By Gender") + theme(plot.title = element_text(hjust = 0.5))
plot5

# figure A3
plot5 <- plot5 + 
  scale_color_manual(values = rep("black", 9))+
  theme(legend.position = "none")
plot5

# analysis 6: vote by gender
amces6 <- cj(d_match_gender, outcome_q1 ~ anticorruption + ideology + party + bribery + embezzlement + competence + age + education, id = ~id, estimate = "amce", by = ~gender)
diff_amces6 <- cj(d_match_gender, outcome_q1 ~ anticorruption + ideology + party + bribery + embezzlement + competence + age + education, id = ~id, estimate = "amce_diff", by = ~gender)
plot6 = plot(rbind(amces6, diff_amces6)) + ggplot2::facet_wrap(~BY, ncol = 3L) + ggtitle("Electoral Support: By Gender") + theme(plot.title = element_text(hjust = 0.5))
plot6

# figure A4
plot6 <- plot6 + 
  scale_color_manual(values = rep("black", 9))+
  theme(legend.position = "none")
plot6

